Dynamics of liquid 4 He in confined geometries from Time-Dependent Density 

Functional calculations 
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We present numerical results obtained from Time-Dependent Density Functional calculations of 
the dynamics of liquid 4 He in different environments characterized by geometrical confinement. The 
time-dependent density profile and velocity field of 4 He are obtained by means of direct numerical 
integration of the non-linear Schrodinger equation associated with a phenomenological energy func- 
tional which describes accurately both the static and dynamic properties of bulk liquid 4 He. Our 
implementation allows for a general solution in 3-D (i.e. no symmetries are assumed in order to 
simplify the calculations) . We apply our method to study the real-time dynamics of pure and alkali- 
doped clusters, of a monolayer film on a weakly attractive surface and a nano-droplet spreading on 
a solid surface. 



PACS numbers: 



.10.-m , 68.45.-v , 68.45.Gd 



I. INTRODUCTION 

Density Functional (DF) methods jjj have become in- 
creasingly popular in recent years as a useful computa- 
tional tool to study the properties of classical and quan- 
tum inhomogeneous fluids, especially for large systems 
where the DF methods provides a good compromise be- 
tween accuracy and computational cost. In particular, 
a quite accurate description of the T = properties of 
liquid 4 He has been obtained within a DF approach by 
using the energy functional proposed in Ref. O and later 
improved in Ref. ||. In general, the Density Functionals 
for 4 He depend on a number of phenomenological param- 
eters which are adjusted to reproduce bulk experimen- 
tal properties of liquid 4 He at saturated vapor pressure. 



The resulting functionals are found to describe accurately 
also the properties of the inhomogeneous systems (e.g. 
the liquid-vapor interface of the free 4 He surface) (for a 
thorough comparison between these and other function- 
als used to describe liquid 4 He at T — 0, see Ref.Jij). 

This phenomenological functional has been widely 
used in a variety of problems involving inhomogeneous 
4 He systems like pure and doped clusters J|, Q], elec- 
tron "bubbles" in bulk 4 He||, alkali atom adsorption on 
the surface of liquid 4 He ||], vortices in 4 He clusters JTo|, 
adsorption of liquid 4 He on solid surfaces [pd| fl2| , etc. 

According to Ref. || the static properties of liquid 4 He 
at T=0 are described by the following energy density 
functional (atomic units will be used throughout our pa- 
per): 



Eo[p] =^m j dr(V^f + ± J dv J dv' p(v)p(v')V He . He {\v - r'|) + j dvE c (r) 



(1) 



The first term on the right hand side is the quantum pres- 
sure (corresponding to the kinetic energy of a Bose gas 
of nonuniform density). The second term contains a two- 
body pair potential Vue-He screened at short distances, 
while the last term accounts for correlation effects due to 
the short-range part of the He-He interaction. The ex- 
plicit form of the functional is given in Ref.j^], to which 
the reader is referred for more details (see also the Ap- 
pendix in the present paper). 

The minimization of the energy functional (Q) with re- 
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spect to density variations, subject to the constraint of 
a constant number of 4 He atoms, J p(r)dv — N, leads to 
the equilibrium profile p(r), thus allowing to study the 
static properties of the 4 He system. 

In order to study the dynamics of 4 He an additional 
term (which has no effect on the static properties derived 
from the minimization of the energy functional (1)) has 
been added || to the functional (1) in the form of a 
phenomenological current term: 



H J = P M M y(rr-^Jdr'V J (\r. 



r'|)p(r)p(r')[v(r)-v(r')] 5 



(2) 
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The first term is the usual hydro-dynamic current den- 
sity, while the second term accounts in a phenomenologi- 
cal way for non-local effects due to the " backflow" current 
density^. The resulting DF (the so-called Orsay-Trento 
Functional), which will be used in our calculations, has 
thus the following form: 



E[p,v] =E [p}+ / dvH, 



(3) 



One appealing feature of this functional, which will 
prove to be essential in the time-dependent calculations 
presented in the following Sections, is that it reproduces 
not only the static structure factor but also the bulk dis- 
persion relations of sound excitations in liquid 4 He, i.e. 
the effective current-current interaction Vj in (2) is em- 
pirically adjusted in order to reproduce the phonon-roton 
dispersion curve of bulk 4 He and the dynamic structure 
function as well. 

The extended functional (3) has been applied in 
the past to the study of density oscillations in 4 He 
clusters [|l3), of surface excitations in 4 He films || || and 
of other related dynamical phenomena |Q . In all of the 
above applications a linear response approximation was 
used, as well as suitable symmetries were assumed in or- 
der to make the calculations feasible. 

There are many different situations, however, where 
a more general approach to investigate the microscopic 
dynamics of liquid 4 He in arbitrary environments, often 
characterized by geometrical confinement, is highly desir- 
able. For instance, modifications of the excitation spec- 
tra due to geometrical confinement are expected in the 
case of 4 He films adsorbed on solid surfaces. A second 
example involves impurity-doped 4 He clusters. Atomic 
impurities are currently used as experimental probes of 
the superconducting behavior of 4 He clusters: the cou- 
pling of the impurity motion with the cluster dynamics is 
crucial for a comprehension of the observed spectra. As a 
third example, we mention the problem of the dynamics 
of wetting phenomena: weakly attractive surfaces like the 
heavy alkali metal surfaces are not wet by liquid 4 He at 
very low temperature, and thus finite amounts of 4 He on 
these surfaces will form a droplet characterized by a non- 
zero contact angle. Even on a highly uniform surface, 
however, the contact angle between the 4 He droplet and 
the surface is extremely hysteretic and its value strongly 
depends on whether the contact line is advancing or re- 
ceding. Informations from microscopic calculations of the 
dynamic behavior of a liquid droplet /film spreading on 
a weak solid surface can represent a valuable help to un- 
derstand these effects. 

We present here calculations, done in the framework 
of Time-Dependent Density Functional (TD-DFT) phe- 
nomenological theory, of the dynamical properties of He 
(at T=0) in confined geometries. Our methods allows to 
calculate the direct real-time evolution of a 4 He system 
in arbitrary 3-D geometries. Some applications of the 
method, showing its efficiency and capabilities, will be 
presented in the following Sections. 



II. TIME-DEPENDENT DENSITY 
FUNCTIONAL THEORY 

In order to extend the DFT scheme described in the 
previous Section to the domain of time evolution (Time 
Dependent Density Functional Theory, TD-DFT), one 
can use a variational principle as follows. One starts 
from a suitable action integral A written in terms of a 
complex "wave function" \f(r, t) 



dtdr <{ £ - «**-^- 



(4) 



where the energy density £ is defined by E\p, v] = J dr£. 
From the condition of stationarity 



a time-dependent Euler-Lagrange equation follows: 
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■U\p,y] 



* = H^> = i- 



.9* 

dt 



(5) 



(6) 



The "effective potential" U appearing in the above non- 
linear time-dependent Schrodinger equation is defined 
as the variational derivative of the energy functional, 
U[p, v] = SE[p(r)]/S^f* and its explicit expression is 
given in the Appendix. From the solution \E' = 0e* e of 
Eq.(^|) one can get the time-dependent density p(r,t) — 
4> 2 and the fluid velocity field v(r,t) = ;jjV9. 

We describe in the following Section the numerical pro- 
cedure that we used in order to integrate numerically the 
Eq.®E|. 



III. COMPUTATIONAL SCHEME 

The numerical solution of eq.(|^), i.e. the wave function 
at an arbitrary time t, ^(r, t), is obtained by using the 
Crank-Nicholson's scheme (CN), which enables to evolve 
<£(£) to a later time t + At: 
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At) ®(t + At) = l~i—H{t) ^(t) 

(7) 

The CN algorithm is one of the many approximate meth- 
ods that can be used to integrate numerically a time- 
dependent Schrodinger equation like the one shown in 
Eq.(^J). It is 2 nd order accurate in At, as can be easily 
seen by expanding the two sides of Eq.(7) in powers of 
At. Among the advantages of such scheme is its numeri- 
cal stability for long-time evolutions, and the unitariness 
(i.e. the number of particles is conserved during the time 
evolution) |^|. The CN recursion can be easily solved 
iteratively for each time step evolution. We find that 
three iterations are usually enough to guarantee a stable 
evolution even during long simulation times. 

The application of the CN formula requires the calcu- 
lation of the action of the " Hamiltonian" H = — V 2 + 
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U[p, v] to the complex wavefunction A direct eval- 
uation of the potential term U (whose explicit expres- 
sion is reported in the Appendix) would be prohibitively 
costly, due to the presence of convolution integrals of the 
form I(r) = J dr'f(r — r')g(r'), whose evaluation in real 
space requires a huge number of floating point opera- 
tions. For such reason we decided to employ a mixed 
real-reciprocal space representation of the basic quanti- 
ties ^ and of the density p. Since the Fourier transform 
of a convolution is the product of the Fourier transforms 
of the two terms entering the integrals, one can perform 
efficiently the calculations of a convolution integral on 
a computer by using Fast Fourier Transforms (FFT). 
The use of FFT implies Periodic Boundary Conditions 
(PBC), and thus our calculations have been performed 
using a periodically repeated supercell containing N He- 
lium atoms. Both \I/ and p are expanded in plane waves: 
\I/(r) = ^ G $Ge !G r , and similarly for p. The G's arc 
the supercell wave-vectors, G = Tr(n/L x ,m/L y ,p/L z ), 
with n,m,p = 0, ±1, ±2, ... and L x , L y , L z are the sides of 
the supercell (a primitive orthorhombic supercell is used 
in all our calculations). The number of plane waves in 
the above expansions is chosen such as to give converged 
values for the total energy E and for the structural pa- 
rameters of the 4 He system under investigation. In the 
various applications that we report in the following, the 
supercell size was always chosen in such a way to make 
negligible the interaction of the 4 He system with its re- 
peated images. 

The action of the kinetic energy operator T = — ^jgV 2 
on the wavefunction is also evaluated in reciprocal 
space, i.e. 

f * = Ym fft{gHg} (8) 

Finally, the real-space expression of Hty , computed as 
described above, is used in the iteration formula of the 
CN scheme to give the evolved wave- function ^(t + At). 
A new potential U is thus computed from the updated 
density and velocity, and the procedure is repeated for 
the next time-step. Typical time steps At which arc 
found to give stable time evolutions are of the order of 
1-3 femtoseconds. 

We stress the fact that no geometric or simmctry con- 
straints are imposed during the time evolution to reduce 
the computational effort. Our computational scheme is 
thus general and can be applied in particular to those 
situations characterized by geometrical confinement of 
4 He. We will give in the following Sections few examples 
of such case. 

IV. 4 HE BULK EXCITATIONS 

One important feature of the most recent implementa- 
tion of DF for 4 He used here is the inclusion of a "back- 
flow" current term (the second term in Equation (2)) ex- 
plicitly depending on the fluid velocity as well as on the 



fluid density, such that the resulting functional is able 
to correctly reproduce the bulk excitation spectrum of 
liquid 4 He. This can be immediately verified within our 
computational scheme as follows. 

We setup an initial state corresponding to a modula- 
tion of the bulk uniform density with a wave vector q: 

<f(r,t = 0) = {p {l + e S m(q - r)}} 1/2 (9) 

where e is a small amplitude and po — 0.02184 A -3 is 
the saturation density at zero temperature and pressure. 
The system is then allowed to evolve in time according 
to the CN scheme described above. By monitoring the 
periodic oscillations of p(r, t) during the simulation, we 
can immediately calculate the mode frequency: 

UJ = Y = w(q) - (10) 

Here T is the period of the observed oscillation. By vary- 
ing the wave vector q we get the whole dispersion curve 
of 4 He, which we compare in Figure 1 with the exper- 
imental data. As it appears from Fig.l, the agreement 
with the measured data is very good in a wide range of 
values of the wavevector q. 

We finally note that if the bulk dispersion relation is 
calculated by neglecting the "backfiow" term in eq.(2), 
results deviating by as much as 50% from the experi- 
mental curva are obtained for q > 0.5 A -1 . 

V. DYNAMICS OF 4 HE CLUSTERS 

The dynamics of liquid 4 He clusters have been studied 
theoretically in recent years by means of a number of dif- 
ferent methods and thus represents a convenient bench- 
mark for our numerical scheme. 

We have considered the monopole and quadrupole vi- 
brational modes of cluster of different sizes, characterized 
by the angular momentum quantum number I = and 
1 = 2, respectively. The dipole mode (1 = 1) corresponds 
to a translation of the cluster as a whole, and should 
thus lie at zero energy. To study the monopole mode, 
which is a radial breathing oscillation of the cluster den- 
sity, we compute the main frequencies by monitoring the 
periodic variations with time of the mean square radius 
of the cluster: 



R(t) = ^Jr*p(r,t)dr. (11) 

To study the quadrupole oscillations we analyze in- 
stead the time evolution of the quadrupole moment ten- 
sor Qif 

Qij(t) = J p(r, t)[3r irj - r 2 S l3 }d 3 r. (12) 

Suitable initial states for both modes are constructed by 
applying to the ground state wave function "to = ^Jp e q 
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(where p eq represents the equilibrium density profile for 
the cluster, calculated as described in Section 1 by direct 
minimization of the energy functional (1)) an excitation 
operator e^ F , where F — r 2 to excite monopole oscilla- 
tions, while F = r 2 y 2 o to excite quadrupole oscillations: 



*(r,t = 0) = e iSir *o(r). 



(13) 



This initial state is thus allowed to evolve according to 
the equation (6). From a frequency analysis of the cal- 
culated M(t) and of Qij(t) we then calculate both the 
amplitudes and the frequencies of the observed oscilla- 
tions. 

We report in Fig. 2 the calculated frequencies of the 
main oscillations observed in the discrete region of the 
spectrum (i.e. below the evaporation threshold ui < |//|, 
H being the chemical potential of a single 4 He atom). 

For comparison, we show in the same Figure the results 
taken from Ref.|l^], obtained within the Random Phase 
Approximation (RPA) and using a zero-range Density 
Functional to describe the 4 He system. The approach 
used in Ref.|ll|, which neglects completely the short- 
range He-He correlation effects, is accurate only when 
applied to the study of excitations involving variations 
in space larger than the average distance between the 
4 He atoms: this is not, however, a major limitation for 
the long-wavelength modes like the ones investigated here 
(external field proportional to r l Yio(r). 

In the same plot we compare also our calculated fre- 
quencies with those obtained from more recent calcu- 
lations done by using a finite-range DF including non- 
local effects, very similar to the one used here, and again 
within RPAQ. The results are not much different from 
those obtained using a zero-range functional(l^], show- 
ing that the effects of non-localities are not important 
in this case due to the q = character of the multipole 
excitations studied here. 

Our calculated values, as it appears from Fig. 2, do 
not differ much from the two previous theoretical re- 
sults, even if they are obtained by explicitly including 
the "backflow" term (see Eq.(2)) in the 4 He density func- 
tional. The reason is that this term becomes important 
in the region where the bulk dispersion relation devi- 
ates from its phonon-like behavior, whereas in the present 
case only low-q excitations are considered. In the next 
sections we will present cases where the inclusion of the 
backflow term is instead crucial in order to get accurate 
results. 



VI. EXCITATIONS OF A 4 HE MONOLAYER 
ADSORBED ON A WEAKLY ATTRACTIVE 
SURFACE 

The study of the structure, growth and excitations of 
liquid 4 He films adsorbed on solid surfaces is an impor- 
tant field of current research on superfluid properties. 



Reliable theories of inhomogeneous liquid 4 He may rep- 
resent a useful complementary tool for studying the prop- 
erties of liquid 4 He interacting with a solid surface [1^] . 

Because of the extremely weak He-He interaction, liq- 
uid 4 He in contact with almost any substrate spreads to 
form a continuous film over the surface, so that vapor and 
substrate are never in contact. A remarkable exception 
to this behavior is found when He is adsorbed on heavy 
alkali metal surfaces: due to the large electron spill-out at 
these surfaces, alkali metals provide the weakest adsorp- 
tion potentials in nature for He atoms. Based on this 
observation, it was suggested |20, Ell that 4 He might 



not wet some heavy alkali metal surfaces. Subsequent 
experiments have confirmed this remarkable prediction 

We study here the excitations of a single monolayer of 
liquid 4 He, adsorbed on a weakly attractive surface. We 
choose here the surface of Rb: this is because Rb is the 
weakest surface that it is actually wet by 4 He at T ~ 0, 
i.e. the stable state of liquid 4 He on Rb is represented 
by a thin liquid film covering the whole surface, rather 
than being represented by a droplet as in the case of 
Cs. The extremely weak He-surface interaction on the 
other hand is expected to act as a small perturbation 
on the film. Thus the excitations studied here should be 
representative to some extent of those for an isolated 4 He 
film. 

To model the interaction with a Rb substrate, we use 
a binding potential V s (z) which describes the interaction 
between Rb, occupying the half space z < 0, and one 4 He 
atom located at a distance z above the surface, which 
is taken to be ideally flat. The explicit form for V s (z), 
originally proposed in Ref . |p7| , has been later revised in 
Rcf . |12] in order to correctly reproduce the experimental 
wetting properties of the 4 He/Rb system. We do not in- 
clude any corrugation of the surface on the atomic scale 
to mimic its actual microscopic structure. This is a good 
approximation for the case of 4 He adsorption on alkali 
metal surfaces, because the experiments |2q| indeed in- 
dicate that a surprisingly smooth surface is seen by ad- 
sorbed 4 He atoms. 

The total energy functional of liquid He interact- 
ing with a Rb surface is thus described by the en- 
ergy functional (3) augmented by the additional term 
J drp(r)V s (r), which describes the interaction energy due 
to the presence of the surface. We have preliminarly stud- 
ied the equilibrium density of liquid 4 He on the Rb sur- 
face, at different areal coverages. Our results are reported 
in Fig. 3, where the density profiles in the direction per- 
pendicular to the surface (which is located at z = 0) are 
shown. Note that as the coverage increases, a second 
layer starts to form above the first layer. Quite arbi- 
trarily, we take as representative of a "monolayer" the 
profile shown in Fig. 3 with a thicker line, corresponding 
to n = 4.2 x 10- 2 A- 2 . 

We thus prepare our monolayer in a suitable non- 
equilibrium initial state, characterized by a density mod- 
ulation with a wave-vector q parallel to the surface. We 
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note that, due to the presence of a symmetry-breaking 
external potential (the He-surface potential) which de- 
pends on z (the coordinate normal to the surface), such 
initial state will be immediately coupled, during its time- 
evolution, to excitations of the film polarized perpendic- 
ularly to the surface plane. 

From a Fourier analysis of the density variations during 
the time evolution from this initial state, we are able to 
get the main excitations of the 4 He film. Fig. 4 shows 
our results. The filled and empty squares show the main 
features observed in the frequency spectrum of the time- 
dependent density p(r, t). The filled squares indicate the 
most intense peaks, whereas the open squares indicate 
the somewhat weaker features in the spectrum. 

For comparison we also show with solid line the ex- 
perimental "ripplon" dispersion relational . It appears 
from Fig. 4 that at low q values the lower disperion curve 
becomes softer that the ripplon curve. This is not sur- 
prising since the ripplon dispersion, which relates to the 
surface excitations of a infinite half-space occupied by 
the liquid, is altered by the finite thickness of the film. A 
simple model calculation 1 30 gives in fact for the oscil- 
lation frequency of a thin ^3e film of thickness a in the 
low-g limit: 



w(q) ~ LOrip P ion{q)(qa) 1/2 



(14) 



Since io r ippion ~ <7 3 ^ 2 , the expected dispersion scales as 
q 2 , in qualitative agreement with our findings. Although 
a definite character cannot be assigned unambiguously 
to the observed mode shown in Fig. 4 due to the strong 
mixing between the two bands, we suggest that the up- 
per branch should have a strong component of density 
modulations perpendicular to the surface. To verify this 
we have induced in the monolayer a simple q = mode 
where the 4 He film is initially rigidly shifted (by a very 
small amount) towards the surface. Due to the presence 
of the holding He-surface potential, the film once let free 
to evolve, execute a motion where it executes oscillations 
mainly perpendicul to the surface, with a main frequency 
at about 5 K: this value is shown with a filled dot at q = 
in Fig. 4. 

We can compare our results for the film excitations 
with those obtained in Ref. (26) by means of a microscopic 
statistical theory, where the dynamic excitations of thick 
4 He films on a Cs surface were considered. In that work 
it is found that ripplon modes appear in the calculated 
spectra, whose dispersion follows the q 3 ^ 2 law for sur- 
face excitations, as well as " interfacial" ripplons local- 
ized mainly at the 4 He/surface interface. In particular, 
the lowest mode calculated in Ref. |^6| (corresponding the 
the lower branch in Fig. 4) has a ripplon character. The 
second mode has a "longitudinal" character (i.e. parallel 
to the surface) on stronger substrates, only on very weak 
surfaces (i.e. Cs) it has a visible perpendicular compo- 
nent. A linear (third sound) dispersion has been also 
observed in Ref. Eg] as long as the wavelength is large 
compared with the film thickness. We do not observe 



any such modes, probably because of the finite thickness 
of our system. 

We believe that the upper branch in Fig. 4 is the 
counterpart of a similar excitations observed in neutron 
diffraction experiments on 4 He adsorbed on a graphite 
substrate|nj, where, among other excitations, almost 
dispersionless modes at low frequencies have been ob- 
served. These mode have been interpreted [ |32[ , by means 
of microscopic calculations, as standing waves of the 4 He 
liquid perpendicular to the substrate. We cannot com- 
pare directly our results with the experimental measure- 
ments of |nj and with the microscopic calculations of 
Ref. Q because of the different substrate used (graphite 
is a much "stronger" substrate than alkali metal surfaces) 
and because multylayer films were considered in these ref- 
erences. We hope that our results will stimulate further 
experimental measurements, in particular to study the 
excitations of 4 He monolayers on weakly attractive sur- 
faces. 



VII. 



DYNAMICS OF ALKALI-DOPED 4 HE 
CLUSTERS 



The spectroscopy of doped liquid 4 He clusters has be- 
come a current tool for the study of superfluidity in 4 He 
droplets As a consequence of the ultra- weak alkali- 
He interaction, an alkali atom picked up by a 4 He cluster 
has its stable state in a "dimple" on the surface of the 
cluster^, rather than being solvated into its inte- 
rior like most of the other atomic impurities. Thus alkali 
could provide an excellent probe for the surface excita- 
tion of superfluid nanodroplets ]34[] . 

The analysis and interpretation of experimental studies 
of the formation of 4 He- alkali complexes in doped clusters 
would benefit from theoretical work of the excitations of 
impurities attached to 4 He clusters. 

Although for light alkali the motion of the impurity 
can be considered to a first approximation as occurring 
in a static 4 He environment (i.e. He does not adjust ap- 
preciably to the istantaneous atom position), however, 
for heavy alkali (Rb,Cs) a full dynamical approach is re- 
quired, where the 4 He atoms in the cluster are allowed 
to adjust dynamically upon the motion of the adsorbed 
impurity. 

We address here the problem of finding the spectrum of 
vibrational excitations of a "dimple" state. We consider 
in particular the system composed of a single Rb atom (in 
its electronic ground-state) attached to a 300-atom 4 He 
cluster. The Rb-He interaction is taken from Ref. ]3l| . 

Fig. 5 shows the calculated equilibrium density profile 
(in a plane containing the center of the cluster and the 
Rb atom) for the " dimple" state representing the lowest 
energy configuration of the Rb- 4 He system. 

By applying to the Rb atom a small initial momen- 
tum towards the surface of the cluster, we observe that 
the impurity starts oscillating around its equilibrium po- 
sition shown in Fig. 5, and at the same time the He 
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density changes in time to adjust dynamically to the is- 
tantaneous Rb position. Our results for the Rb dynamics 
are shown in Fig. 6. In the upper panel we show the cal- 
culated istantaneous position of the Rb atom as a func- 
tion of time (measured with respect to the Rb-cluster 
center-of-mass). A straightforward Fourier transform of 
the signal shown in the upper panel allows to compute 
the frequency spectrum shown in the lower panel. 

We compare our results with a previous estimate Q 
based on static DF calculations, where the frequency of 
the oscillations of a Rb atom attached to a 300-atom 
cluster was estimated approximately by assuming that 
the He atoms adiabatically adjusted to the istantaneous 
position of the impurity. In that case the value of the fre- 
quency was extracted from the shape of the impurity- He 
potential energy surface, giving uj ~ 0.7 K. It appears 
from our results that a more complex dynamics results, 
due to the coupling of the Rb atom motion with the sur- 
face oscillations of the cluster. The frequency of a sur- 
face m ode of an N-at om 4 He cluster can be written as 
lj = y/l(l- l)(l + 2)u /VN where lu = 3A5K. This 
is approximately ~ 0.6 K for the lowest I = 2 mode, i.e. 
close to the natural Rb frequency calculated in Ref. |36j . 
A large amount of mixing is thus expected, which deter- 
mines the appearance of the complex vibrational spec- 
trum shown in Fig. 6. 

Experimental measurements of such oscillations are 
under way |37]] to confirm our predictions. A more com- 
plex and challenging application of time-dependent meth- 
ods, i.e. the study of the oscillations of an alkali impurity 
in an excited electronic state, is currently in progress. 



VIII. SPREADING OF A 4 HE NANODROPLET 
ON A WEAKLY ATTRACTIVE SUBSTRATE 

The mechanism underlying the spreading of a liquid 
droplet on a solid substrate is a long-standing problem. 
Experiments done on a number of differente (classical) 
molecular fluids have revealed universal laws such as 
the Tanner' spreading law [ |38| , where the droplet radius 
grows as r(t) oc t 1 / 10 , as well as fascinating phenomena 
on the atomistic scale (such as the appearence of a mono- 
layer film advancing in front of the dropj39|, two differ- 
ent regimes of the growth of the radius with timejfO), 
dynamical "layering", etc.). Other experiments ]39|] have 
shown that one or more monomolecular precursor layers 
spread ahead of the droplet cap with an average radius 

R(t) ~ y/t. 

Numerical studies of droplet spreading have been per- 
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formed 1 40 , [41|, |42|, |43| in recent years. Despite many 
efforts, however, many microscopic details of the spread- 
ing process are however still to be known. 

Even more challenging is the behavior of a superfluid 
4 He droplet on a weakly attractive surface: it has been 
found that 4 He droplets adsorbed on a Cs substrate have 
spreading and flow properties that are not simple conse- 
quences of bulk superfluid behavior pi. When a 4 He 



droplet is deposited on a Cs surface (the only mate- 
rial known that is not wetted by superfluid 4 He at very 
low temperature), even on a highly uniform surface the 
contact angle between the droplet and the surface is 
extremely hysteretic and its value strongly depends on 
whether the contact line is advancing or receding. Su- 
perfluid 4 He droplets on Cs are also remarkable because 
they can resist flow against a substantial chemical poten- 
tial gradient Q. 

We have studied the spreading of a 200-atom He cluster 
on a weakly attractive surface. We choose for the latter 
the case of Cs (but very similar results are obtained for a 
Rb surface). To simulate the Cs surface, we use the "ab 
initio" potential developed in Ref. [^7| , which has applied 
successfully to the wetting properties of 4 He an a Cs sur- 
face |l2]. The initial state is represented by a spherical 
cluster (whose density profile is obtained by minimization 
of the static functional (2) in the absence of the attrac- 
tive He-surface potential), placed in the vicinity of the 
surface. The downward force responsible for the spread- 
ing of the droplet is provided solely by the long-range 
Van der Waals forces exerted by the substrate (due to 
the microscopic size of the droplet, the effect of gravity 
on spreading is negligible). 

The spreading will continue, until the stable state is 
eventually reached. The nature of the latter depends on 
the strength of the He-surface potential: in the case of 
Cs, which is not wet by He at T=0, the final configuration 
will be that of an almost spherical cap characterized by a 
contact angle of ~ 40 ° (partial wetting behavior) Q . 
More attractive substrate, on the other hand, will eventu- 
ally be covered by a uniform microscopic film (complete 
wetting behavior) . Due to the finiteness of our systems 
and our use of periodic boundary conditions our simu- 
lations are necessarily interrupted when the droplet hits 
its repeated images and starts coalescing with them. 

The sequence characterized the spreading process is 
shown in Fig. 7. The overall duration of the simulation 
shown in the Figure is about 2ns. Interestingly enough, 
a precursor layer seems to form under the cap of the 
spreading cluster. We have computed from the cluster 
density profile p(r, t) the average radius R of this precur- 
sor layer. The dependence of R from the time, shown in 
Fig. 8, shows a "fast" regime where the precursor layer 
expands linearly in time, as found in similar simulations 
for classical fluids^]. Then a slowing down occurs, in 
correspondence of the formation of the second layer (see 
panel (f) in Fig. 7). It seems that after a transient, an 
almost linear spreading occurs again, although with a 
somewhat lower velocity. The speed from the linear por- 
tion of Fig. 8 gives a velocity of about ~ 50 m/s, which is 
close to the critical Landau velocity. The slowing down 
is perhaps due to the spontaneous formation of vorti- 
cal rings during the spreading. One such ring is shown 
in Fig. 9 by means of a vector plot showing the current 
density J = p(r)v(r) in a plane containing the center of 
the cluster and perpendicular to the Cs surface. Fig. 9 
also clearly shows the flow of atoms from the droplet cap 
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feeding the underlying precursor film. 

Although transient precursor films spreading linearly 
with time has been obtained from MD simulations of 
classical fluids ^3|, they have not been seen in any 
experiment. This is probably due to the fact that they 
represent the very early stage of the spreading process: 
as the equilibrium shape of the droplet is reached, the 
initial fast spreading is found to slow down to a power 
law R(t) ~ t 1 / 7 ]^, when a regime dominated by dissipa- 
tion due to the viscous flow starts to dominate. It would 
be extremely interesting to study how this slowing down 
(which is only suggested from our short-time results in 
Fig. 8) occurs in the case of superfluid flow. 

The presence of surface disorder is expected to alter 
the dynamics of liquid 4 He films spreading on a solid sur- 
face, because the three-phase contact line of the spread- 
ing film/droplet can easily be pinned to these defects, 
modifying the velocity of the contact line as it advances 
(or recedes) on the surface. We are currently investi- 
gating the effect of simple defects on the atomic-scale 
(adatoms or vacancies) on the droplet spreading process 
discussed above. 



geometrical confinement. We have shown the feasibility 
and efficiency of fully 3-D TD-DFT calculations to study 
the dynamics of liquid He-4 in a number of representa- 
tive situations: the dynamics of pure and alkali-doped 
4 He clusters, the excitations of a liquid 4 He monolayer 
adsorbed on a weakly attractive surface and the dynam- 
ics of a 4 He nano-droplet spreading on a Cs surface. 
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IX. SUMMARY 

We have presented TD-DFT calculations of the dy- 
namical properties of 4 He in situations characterized by 



Appendix 

The effective potential entering Eq.(6) can be read- 
ily evaluated by functional differentiation of the energy 
functional (3) (see Ref.pfl for its detailed expression): 



U[p,v] = V s (r)+ I dr , p(v')Ve(\r-r'\) + ^p(r) 2 + -^p(rf+ I dr' p(r')[c 2 U h (\r - r'|)p(r') + C3 n h (|r - r'|)p 2 (r')] 

+ Sf i 1 ~ / dT ' i 1 ~ Vr ' P(r,) ' VrF(|r ' ~ r|) " Y / dr ' Vji]r ~ AW,t){v{r) - v(v')f 

+ ^)V-|dr'g(r,r') (15) 



V s (r) represents an additional external potential act- 
ing on the 4 He system. The second term contains a two- 
body 4 He- 4 He pair potential Ve(r) screened at distances 
shorter than a characteristic length he, while the third 
and the fourth term (correlation terms), which contains 
the average of the density over a sphere of radius he, p r , 
accounts for the internal kinetic energy and for the in- 
creasing contribution of the hard-core He-He repulsion 
when the density is increased. The last two terms rep- 
resent the contribution to the potential U of the "back- 
flow" term (2) introduced in Ref. || in order to reproduce 
the bulk excitation spectrum of liquid 4 He. In particu- 



lar, the last term contains the "backflow" current Js = 
/ dr'g(r, r') = / dr'Vj(\r - r'|)p(r)p(r')(v(r) - v(r')) 

The free parameters he, C2, C3, a s are adjusted in order 
to reproduce a number of experimental properties of bulk 
liquid 4 He. For a detailed description of the various terms 
and the numerical values of the adjustable parameters, 
we refer the interested reader to Ref.pl. 
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Figure Captions 

Figure 1: Bulk 4 He excitation spectrum. The points 
show the calculated spectrum, while the solid line shows 
the experimental dispersion relation. 

Figure 2: Calculated monopolc (upper lines) and 
quadrupole (lower lines) oscillation frequencies of liquid 
4 He clusters, shown as a function of the cluster sizes. 
Dots: this work; squares: results from Ref. [jDl ; triangles: 
results from Ref. || . 

Figure 3: Equilibrium density profiles for 4 He films 
adsorbed on a Rb surface. The z direction is perpen- 
dicular to the surface plane. Different curves represent 
increasing values of the surface coverages (number of 4 He 
atoms per unit surface area), from n = 1.6 x 10~ 2 A~ 2 
(lower profile) to n = 7.2 x 10~ 2 A~ 2 (upper profile). 
The thicker line (defining our "monolayer") is drawn at 
n = 4.2 x 1CT 2 A- 2 . 

Figure 4: Calculated dispersion curves for a 4 He mono- 
layer adsorbed on a Rb surface (squares) . The dotted line 
is the bulk dispersion relation, while the solid line shows 
the experimental "ripplon" dispersion. 

Figure 5: Density contour plot for a 300-atom 4 He 
cluster with a Rb atom attached on its surface: the dot 
shows the equilibrium position of the Rb impurity. 
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Figure 6: Upper panel: Rb displacement with respect 
to the center-of mass of the Rb-cluster system, measured 
along the radial direction, as a function of the simulation 
time. Lower panel: Fourier spectrum (in arbitrary units) 
of the time series shown in the upper panel. 

Figure 7: Sequence of 4 He density contour plots show- 
ing the spreading of a 200-atom cluster on a Cs surface. 



Figure 8: Radius of the precursor layer as a function 
of time. 



Figure 9: Vector plot of the current density pv in a 
plane perpendicular to the surface and passing through 
the center of the cluster. The configuration shown corre- 
sponds to the density contours in panel (e) of Fig. 7. 
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